# Program used for Monte-Carlo results proposed in the paper #
##############################################################

rm(list=ls(all=TRUE))


# Procedure used to get results as tex table

source("C:\\print_results.R")


# Library containing the optimization algorithm

library(limSolve)

alphaC	<- 0		 		# Constant
alphaT	<- .3				# Treatment effect



# Data generating process #
###########################

dgpinter <- function(N,N0,T,T0,alphaT,alphaX,alphaC,supp,fei0,fet,vecL0,vecF,link,nbf,rho,Vr,dgp) {

if (dgp<3) {
	if (link=="normal")  { u <- rnorm(N*T,0,1) }
	if (link=="uniform") { u <- (runif(N*T,0,1)-.5)/.5*sqrt(3) }
	if (rho>0) {
		for (i in 1:N) {
		for (t in 2:T) {
			u[(i-1)*T+t]=rho*u[(i-1)*T+t-1]+u[(i-1)*T+t]
}
}
}
}

if (dgp==3) {
if (link=="normal")  { 
		a=(Vr-1)/(Vr+1)
		low <- (runif(N*T,1)<.5)
		u <- (low*sqrt(1+a)+(1-low)*sqrt(1-a))*rnorm(N*T,0,1)
}
}

vecD		<- c(rep(1,N0),rep(0,N-N0))
vecT		<- c(rep(0,T0-1),rep(1,T-T0+1))

fei		<- as.matrix(fei0,N,1)
fei[1:N0,1]	<- fei[1:N0,1] + matrix(supp,N0,1)
vecL		<- as.matrix(vecL0[,1:nbf],N,nbf) 
vecL[1:N0,]	<- vecL[1:N0,] + matrix(supp,N0,nbf)

alphaTb 	<- matrix(alphaT,N*T,1)	
X	    	<- as.matrix(vecD %x% vecT)

y	 	<- X * alphaTb + 
		   alphaC +
		   as.matrix(fei %x% rep(1,T),N*T,1) +
		   as.matrix(rep(1,N) %x% fet,N*T,1) + 
 		   matrix(vecF[,1:nbf] %*% t(vecL),N*T,1) + 
		   u

var <- list(y,X)
names(var)  <- c("y","X")
var
}



# Procedure used to find synthetic control weights #
####################################################

wipop <- function(Z0,Z1) {

E <- matrix(1,1,ncol(Z0))
F <- 1

b <- 1

G <- diag(ncol(Z0))
H <- matrix(0,ncol(Z0),1)

res <- lsei(A=Z0, B=Z1, E=E, F=F, G=G, H=H, type=2)
w <- res$X

w
}



# Estimation of parameters using synthetic controls #
#####################################################

coefsSC <- function(y,X) {

yy0  <- matrix(y[(N0*T+1):(N*T),1],T,N-N0)
y0   <- yy0[1:(T0-1),]

yNT  <- matrix(0,N0*T,1)

for (j in 1:N0) {
y1 <- matrix(y[((j-1)*T+1):((j-1)*T+T0-1),1],T0-1,1)
ww <- wipop(y0,y1)
yNT[((j-1)*T+1):(j*T),1] <- yy0%*%ww
}

w0 <- rep(1,N0)%x%c(rep(0,T0-1),rep(1,T-T0+1))

ATTSC <- t(w0)%*%(y[1:(N0*T),1]-yNT)/(N0*(T-T0+1))

ATTSC
}



# Estimation algorithm steps for least-squares estimations #
############################################################

# Step for parameters of explanatory variables

coefX <- function(yp,Xp,iXpXp,betaL,betaF) {
betaX <- iXpXp %*% crossprod(Xp,yp-matrix(betaF %*% t(betaL),nrow(yp),1))
betaX
}

# Step for constant, individual and time additive effects
# where there is no explanatory variable

coefCit0 <- function(ym,yi,yt) {
coefC <- as.numeric(ym)
coefi <- yi - coefC
coeft <- yt - coefC
coeff <- list(coefC,coefi,coeft)
names(coeff) <- c("coefC","coefi","coeft")
coeff
}

# Step for time interactive effects
# where there is no explanatory variable

coefF0 <- function(yp,nbf) {
ww <- matrix(0,T,T)
for (i in 1:(nrow(yp)/T)) {
wi <- yp[((i-1)*T+1):(i*T)]
ww <- ww+crossprod(t(wi))
}
oww <- order(eigen(ww)$values,decreasing=TRUE)
oww <- oww[1:nbf]
betaF <- matrix(sqrt(T)*eigen(ww)$vectors[,t(oww)],T,nbf)
if (betaF[1]<0) {betaF=-betaF}
betaF
}

# Step for individual interactive effects
# where there is no explanatory variable

coefL0 <- function(yp,betaF) {
betaL <- matrix(0,N,ncol(betaF))
for (i in 1:N) {
betaL[i,] <- crossprod(yp[((i-1)*T+1):(i*T)],betaF)/T
}
betaL
}

# Step for constant, individual and time additive effects

coefCit <- function(ym,Xm,yi,Xi,yt,Xt,betaX) {
coefC <- as.numeric(ym - Xm%*%betaX)
coefi <- yi - Xi%*%betaX - coefC
coeft <- yt - Xt%*%betaX - coefC
coeff <- list(coefC,coefi,coeft)
names(coeff) <- c("coefC","coefi","coeft")
coeff
}

# Step for time interactive effects

coefF <- function(yp,Xp,betaX,nbf) {
N <- nrow(yp)/T
ww <- matrix(0,T,T)
for (i in 1:N) {
wi <- yp[((i-1)*T+1):(i*T)]-Xp[((i-1)*T+1):(i*T),]%*%betaX
ww <- ww+crossprod(t(wi))
}
oww <- order(eigen(ww)$values,decreasing=TRUE)
oww <- oww[1:nbf]
betaF <- matrix(sqrt(T)*eigen(ww)$vectors[,t(oww)],T,nbf)
if (betaF[1]<0) {betaF=-betaF}
betaF
}

# Step for individual interactive effects

coefL <- function(yp,Xp,betaX,betaF) {
N <- nrow(yp)/T
betaL <- matrix(0,N,ncol(betaF))
for (i in 1:N) {
betaL[i,] <- crossprod(yp[((i-1)*T+1):(i*T)]-
             Xp[((i-1)*T+1):(i*T),]%*%betaX,betaF)/T
}
betaL
}


# Estimation of parameters using Bai (2009) #
#############################################

coefsbai <- function(y,X,nbf,initL,initF,iter,expl) {

# Quantities computed once and for all

N <- nrow(initL)
T <- nrow(initF)

numi <- seq(1:N) %x% rep(1,T)
numt <- rep(1,N) %x% seq(1:T)
ym   <- apply(y,2,mean)
yi   <- as.matrix(tapply(y,numi,mean),N,1)
yt   <- as.matrix(tapply(y,numt,mean),T,1)
yp   <- y-as.matrix(yi %x% rep(1,T))-as.matrix(rep(1,N) %x% yt)+ym

if (expl==1) {
Xm   <- t(apply(X,2,mean))
Xi   <- matrix(0,N,ncol(X))
for (k in 1:ncol(X)) { Xi[,k] <- tapply(X[,k],numi,mean) }
Xt   <- matrix(0,T,ncol(X))
for (k in 1:ncol(X)) { Xt[,k] <- tapply(X[,k],numt,mean) }
Xp   <- X-Xi%x%rep(1,T)-rep(1,N)%x%Xt+Xm%x%rep(1,N*T)
iXpXp <- solve(crossprod(Xp))
}

if (expl==0) {
cf     <- coefCit0(ym,yi,yt)
coefF  <- coefF0(yp,nbf)
coefL  <- coefL0(yp,coefF)
}

if (expl==1) {
for (i in 1:iter) {
if (i==1) { coefX <- coefX(yp,Xp,iXpXp,initL,initF) }
else	    { coefX <- coefX(yp,Xp,iXpXp,coefL,coefF) }
cf     <- coefCit(ym,Xm,yi,Xi,yt,Xt,coefX)
coefF  <- coefF(yp,Xp,coefX,nbf)
coefL  <- coefL(yp,Xp,coefX,coefF)
}
}

if (expl==0) {
coeff <- list(cf$coefC,cf$coefi,cf$coeft,coefF,coefL)
names(coeff) <- c("coefC","coefi","coeft","coefF","coefL")
}

if (expl==1) {
coeff <- list(coefX,cf$coefC,cf$coefi,cf$coeft,coefF,coefL)
names(coeff) <- c("coefX","coefC","coefi","coeft","coefF","coefL")
}

coeff
}



# Estimation of standard errors for Bai (2009) #
################################################

stdbai_iid <- function(y,X,coefX,coefC,coefi,coeft,coefL,coefF) {

# Quantities computed once and for all

numi <- seq(1:N) %x% rep(1,T)
numt <- rep(1,N) %x% seq(1:T)
ym   <- apply(y,2,mean)
Xm   <- t(apply(X,2,mean))
ymi  <- as.matrix(tapply(y,numi,mean),N,1)
Xmi  <- matrix(0,N,ncol(X))
for (k in 1:ncol(X)) { Xmi[,k] <- tapply(X[,k],numi,mean) }
ymt  <- as.matrix(tapply(y,numt,mean),T,1)
Xmt  <- matrix(0,T,ncol(X))
for (k in 1:ncol(X)) { Xmt[,k] <- tapply(X[,k],numt,mean) }

yp <- y-as.matrix(ymi %x% rep(1,T))-as.matrix(rep(1,N) %x% ymt)+ym
Xp <- X-Xmi%x%rep(1,T)-rep(1,N)%x%Xmt+Xm%x%rep(1,N*T)

MF <- diag(T)-coefF%*%solve(crossprod(coefF))%*%t(coefF)

llm1 <- solve(crossprod(coefL)/N)

eps  <- matrix(0,N*T,1)
som  <- 0

for (i in 1:N) {
somi <- 0
Xi  <- Xp[((i-1)*T+1):(i*T),]
li  <- matrix(coefL[i,],ncol(coefL),1)

eps[((i-1)*T+1):(i*T),1] <- y[((i-1)*T+1):(i*T),1]-X[((i-1)*T+1):(i*T),]%*%coefX-
				    coefF%*%li-coefC-coefi[i,1]-coeft

for (k in 1:N) {
	Xk  <- Xp[((k-1)*T+1):(k*T),]
	lk <- matrix(coefL[k,],ncol(coefL),1)	
	aik <- as.numeric(t(li)%*%llm1%*%lk)
	somi <- somi + aik*MF%*%Xk/N
}
zi <- MF%*%Xi-somi
som <- som + crossprod(zi)/(N*T)
}
sig2 <- as.numeric(crossprod(eps)/(N*T))

stdbai <- matrix(sqrt(diag(sig2*solve(som))/(N*T)),ncol(som),1)
stdbai
}





# Estimation of parameters using Bai (2009), ATT #
##################################################

coefsbaiATT <- function(y,X,nbf,initL,initF,iter,iterEM) {

# Initialisation of parameters

yy <- as.matrix(y[(N0*T+1):(N*T),1])

cf <- coefsbai(yy,XX,nbf,matrix(initL[(N0+1):N,],N-N0,ncol(initF)),
               initF,iter,0)

ccF <- solve(crossprod(cbind(matrix(1,T0-1,1),cf$coefF[1:(T0-1),])))%*%
                       t(cbind(matrix(1,T0-1,1),cf$coefF[1:(T0-1),]))

coeff <- matrix(0,N0,nbf+1)

for (i in 1:N0) {

coeff[i,] <- t(ccF%*%(y[((i-1)*T+1):((i-1)*T+T0-1),1]-cf$coeft[1:(T0-1),1]-
             cf$coefC))

}

yNT <- matrix(0,N0*T,1)

for (k in 1:iterEM) {

# E-step of EM algorithm

for (i in 1:N0) {

yNT[((i-1)*T+1):((i-1)*T+T0-1),1] <- y[((i-1)*T+1):((i-1)*T+T0-1),1]

if (k==1) {
yNT[((i-1)*T+T0):(i*T),1] <- cf$coefC+cf$coeft[T0:T,1]+
      cbind(matrix(1,T-T0+1,1),cf$coefF[T0:T,])%*%matrix(coeff[i,],nbf+1,1)
}
else {
yNT[((i-1)*T+T0):(i*T),1] <- cf$coefC+cf$coeft[T0:T,1]+
      cbind(matrix(1,T-T0+1,1),cf$coefF[T0:T,])%*%
      matrix(cbind(cf$coefi[i,1],t(cf$coefL[i,])),nbf+1,1)
}

}

yb <- rbind(yNT,yy)

# M-step of EM algorithm

cf <- coefsbai(yb,X,nbf,initL,initF,iter,0)

}

yNT <- cf$coefC+rep(1,N0)%x%cf$coeft+
       matrix(cbind(matrix(1,T,1),cf$coefF)%*%
       t(cbind(cf$coefi[1:N0,1],(cf$coefL[1:N0,]))),N0*T,1)


# Estimation of treatment effect without constraint

numiT <- rep(1,N0) %x% c(rep(0,T0-1),rep(1,T-T0+1))

ATT <- matrix(numiT,1,N0*T)%*%(y[1:(N0*T),1]-yNT)/(N0*(T-T0+1))

# Estimation of treatment effect with constraint

yy0  <- matrix(y[(N0*T+1):(N*T),1],T,N-N0)
y0   <- yy0[1:(T0-1),]

cT   <- cbind(rep(1,T0-1),cf$coefF[1:(T0-1),])
cNT  <- cbind(rep(1,T),cf$coefF)

www  <- matrix(0,N0,N-N0) 

somXX <- matrix(0,(N-N0)*(nbf+1),(N-N0)*(nbf+1))
somXY <- matrix(0,(N-N0)*(nbf+1),1)

for (j in 1:N0) {

y1 <- matrix(y[((j-1)*T+1):((j-1)*T+T0-1),1],T0-1,1)

ww <- wipop(y0,y1)

www[j,] <- t(ww)

XT <- t(ww)%x%cT

somXX <- somXX+crossprod(XT)
somXY <- somXY+crossprod(XT,y1-cf$coefC-cf$coeft[1:(T0-1),1])
}

somXX <- somXX + diag(N-N0)%x%crossprod(cNT)
somXY <- somXY + (diag(N-N0)%x%t(cNT))%*%
                 (y[(N0*T+1):(N*T),1]-cf$coefC-rep(1,N-N0)%x%cf$coeft)

varcoef <- solve(somXX)%*%somXY

varcoef <- t(matrix(varcoef,nbf+1,N-N0))
coefiNT <- matrix(varcoef[,1],N-N0,1)
coefiT  <- www%*%coefiNT
coefi   <- rbind(coefiT,coefiNT)
coefLNT <- matrix(varcoef[,2:(nbf+1)],N-N0,nbf)              
yy      <- y-coefi%x%rep(1,T)-cf$coefC

coefLT <- www%*%coefLNT
coefL <- rbind(coefLT,coefLNT)

# Last step to improve efficiency

yyT    <- matrix(yy[1:(N0*T),],T,N0)
yyNT   <- matrix(yy[(N0*T+1):(N*T),],T,N-N0)

coefiLT  <- cbind(rep(1,N0),coefLT)
coefiLNT <- cbind(rep(1,N-N0),coefLNT)

somXX1 <- solve(crossprod(coefiLT)+crossprod(coefiLNT))
somXX2 <- solve(crossprod(coefiLNT))

F1 <- (yyT[1:(T0-1),]%*%coefiLT+yyNT[1:(T0-1),]%*%coefiLNT) %*% somXX1  
F2 <- yyNT[T0:T,] %*% coefiLNT %*% somXX2

coefF <- rbind(F1,F2)

yNT <- cbind(matrix(1,T,1),coefF)%*%t(cbind(coefiT,coefiLT))
yNT <- cf$coefC+matrix(yNT,N0*T,1)

ATTc <- matrix(numiT,1,N0*T)%*%(y[1:(N0*T),1]-yNT)/(N0*(T-T0+1))

coeff <- list(ATT,ATTc)
names(coeff) <- c("ATT","ATTc")

coeff

}



# Normal kernel #
#################

K     <- function(u,h) {
pi 	<- 3.14159 
Ku    <- exp(-(u/h)*(u/h)/2)/sqrt(2*pi)/h
Ku
}



# Estimation of parameters using Bai (2009), matching #
#######################################################

coefsmatch <- function(y,X,nbf,initL,initF,iter,iterEM) {

# Initialisation of parameters

yy <- as.matrix(y[(N0*T+1):(N*T),1])

cf <- coefsbai(yy,XX,nbf,matrix(initL[(N0+1):N,],N-N0,ncol(initF)),
               initF,iter,0)

ccF <- solve(crossprod(cbind(matrix(1,T0-1,1),cf$coefF[1:(T0-1),])))%*%
                       t(cbind(matrix(1,T0-1,1),cf$coefF[1:(T0-1),]))

coeff <- matrix(0,N0,nbf+1)

for (i in 1:N0) {

coeff[i,] <- t(ccF%*%(y[((i-1)*T+1):((i-1)*T+T0-1),1]-cf$coeft[1:(T0-1),1]-
             cf$coefC))

}

yNT <- matrix(0,N0*T,1)

for (k in 1:iterEM) {

# E-step of EM de algorithm

for (i in 1:N0) {

yNT[((i-1)*T+1):((i-1)*T+T0-1),1] <- y[((i-1)*T+1):((i-1)*T+T0-1),1]

if (k==1) {
yNT[((i-1)*T+T0):(i*T),1] <- cf$coefC+cf$coeft[T0:T,1]+
      cbind(matrix(1,T-T0+1,1),cf$coefF[T0:T,])%*%matrix(coeff[i,],nbf+1,1)
}
else {
yNT[((i-1)*T+T0):(i*T),1] <- cf$coefC+cf$coeft[T0:T,1]+
      cbind(matrix(1,T-T0+1,1),cf$coefF[T0:T,])%*%
      matrix(cbind(cf$coefi[i,1],t(cf$coefL[i,])),nbf+1,1)
}

}

yb <- rbind(yNT,yy)

# M-step of EM de algorithm

cf <- coefsbai(yb,X,nbf,initL,initF,iter,0)

}
cL	<- cbind(cf$coefi,cf$coefL)
dumt	<- matrix(c(rep(1,N0),rep(0,N-N0)),N,1)


# Score

myprobit <- glm(dumt~cL,family=binomial(link="logit"))
score    <- as.matrix(predict(myprobit,type="response"),N,1)

hh    <- (nrow(score))^(-1/5)*1.06*as.numeric(sd(score))

ypred <- matrix(0,N0,T)

for (i in 1:N0) {
ypred[i,] <- matrix(colMeans(
		 matrix(K(matrix(score[i,1],N-N0,1)-score[(N0+1):N,1],hh),N-N0,T)*
		 t(matrix(yy,T,N-N0))),1,T)/
	       mean(K(matrix(score[i,1],N-N0,1)-score[(N0+1):N,1],hh))
}

ypred <- matrix(t(ypred),N0*T,1)

numiT <- rep(1,N0) %x% c(rep(0,T0-1),rep(1,T-T0+1))

matchATT <- matrix(numiT,1,N0*T)%*%(y[1:(N0*T),1]-ypred)/(N0*(T-T0+1))

coeff <- list(matchATT)
names(coeff) <- c("matchATT")

coeff

}



# Estimation of parameters using FGLS first diff #
##################################################

coefs1difFGLS <- function(y,X) {
XX <- cbind(X,rep(1,N)%x%rbind(t(rep(0,T-1)),diag(T-1)))
dX <- matrix(0,N*(T-1),ncol(X)+T-1)
dy <- matrix(0,N*(T-1),1)
for (i in 1:N) {
for (t in 1:(T-1)) {
dy[(i-1)*(T-1)+t,1] <- y[(i-1)*T+t+1]-y[(i-1)*T+t]
dX[(i-1)*(T-1)+t,]  <- XX[(i-1)*T+t+1,]-XX[(i-1)*T+t,]
}
}

idXdX <- solve(crossprod(dX))

D  <- 2*diag(T-1)
for (t in 1:(T-2)) {
D[t,t+1] <- -1
D[t+1,t] <- -1
}

coefX <- solve(crossprod(dX))%*%crossprod(dX,dy)

som <- matrix(0,ncol(X)+T-1,ncol(X)+T-1)
for (i in 1:N) {
som <- som+t(dX[((i-1)*(T-1)+1):(i*(T-1)),])%*%D%*%dX[((i-1)*(T-1)+1):(i*(T-1)),]
}

sig2 	<- as.numeric(crossprod(dy-dX%*%coefX)/(2*N*T-sum(diag(idXdX%*%som))))
Vm1	<- solve(sig2*D)

som1 <- matrix(0,ncol(X)+T-1,ncol(X)+T-1)
som2 <- matrix(0,ncol(X)+T-1,1)

for (i in 1:N) {
som1 <- som1+t(dX[((i-1)*(T-1)+1):(i*(T-1)),])%*%Vm1%*%dX[((i-1)*(T-1)+1):(i*(T-1)),]
som2 <- som2+t(dX[((i-1)*(T-1)+1):(i*(T-1)),])%*%Vm1%*%dy[((i-1)*(T-1)+1):(i*(T-1)),]
}

coeff <- solve(som1)%*%som2
coeff <- coeff[1:ncol(X),1] 
coeff
}



# Procedure to compute simulation results #
###########################################

simul2 <- function(N,N0,T,T0,alphaT,alphaX,alphaC,supp,
			 fei0,fet,vecL0,vecF,iter,iterEM,link,nbf,S,rho,Vr,dgp) {

cbai     <- matrix(0,S,1)
cbaiATT  <- matrix(0,S,1)
cbaiATTc <- matrix(0,S,1)
cATTCSC  <- matrix(0,S,1)
cmatchATT <- matrix(0,S,1)
cdiff    <- matrix(0,S,1)

for (s in 1:S) {

# Data generation

var <- dgpinter(N,N0,T,T0,alphaT,alphaX,alphaC,supp,
                fei0,fet,vecL0,vecF,link,nbf,rho,Vr,dgp)

# Initial value of parameters for iteration algorithm

initL <- matrix(0,N,nbf)

initF <- matrix(0,T,nbf)
for (k in 1:nbf) {
initF[k,k] <- sqrt(T)
}

beta <- coefsbai(var$y,var$X,nbf,initL,initF,iter,1)

cbai[s,] <- beta$coefX

beta <- coefsbaiATT(var$y,0,nbf,initL,initF,iter,iterEM)

cbaiATT[s,]  <- beta$ATT
cbaiATTc[s,] <- beta$ATTc

cATTCSC[s,] <- coefsSC(var$y,0)

betam <- coefsmatch(var$y,0,nbf,initL,initF,iter,iterEM)
cmatchATT[s,] <- betam$matchATT

beta <- coefs1difFGLS(var$y,var$X) 
cdiff[s,] <- t(beta)
}

# TABLE CONTAINING RESULTS

results <- matrix(0,18,1)

# Results for Bai (2009), ATT

results[1,1] <- round(mean(cbaiATT),3)-alphaT
results[2,1] <- round(median(cbaiATT),3)-alphaT
results[3,1] <- round(sd(cbaiATT),3)

# Results for Bai (2009)

results[4,1] <- round(mean(cbai[,1]),3)-alphaT
results[5,1] <- round(median(cbai[,1]),3)-alphaT
results[6,1] <- round(sd(cbai[,1]),3)

# Results for Bai (2009), matching

results[7,1] <- round(mean(cmatchATT),3)-alphaT
results[8,1] <- round(median(cmatchATT),3)-alphaT
results[9,1] <- round(sd(cmatchATT),3)

# Results for Bai (2009), ATT constrained

results[10,1] <- round(mean(cbaiATTc),3)-alphaT
results[11,1] <- round(median(cbaiATTc),3)-alphaT
results[12,1] <- round(sd(cbaiATTc),3)

# Results for synthetic controls

results[13,1] <- round(mean(cATTCSC),3)-alphaT
results[14,1] <- round(median(cATTCSC),3)-alphaT
results[15,1] <- round(sd(cATTCSC),3)

# Results for FGLS 1st difference

results[16,1] <- round(mean(cdiff[,1]),3)-alphaT
results[17,1] <- round(median(cdiff[,1]),3)-alphaT
results[18,1] <- round(sd(cdiff[,1]),3)

results
}



# For Bai (2009), computation of finite sample and asymptotic variances #
#########################################################################

simul3 <- function(N,N0,T,T0,alphaT,alphaX,alphaC,supp,
			 fei0,fet,vecL0,vecF,iter,iterEM,link,nbf,S) {

coeff    <- matrix(0,S,1)	
cbai     <- matrix(0,S,1)
vbai     <- matrix(0,S,1)

for (s in 1:S) {
var <- dgpinter(N,N0,T,T0,alphaT,alphaX,alphaC,supp,
                fei0,fet,vecL0,vecF,link,nbf,0,0,1)

# Valeur initiale des param�tres � estimer

initL <- matrix(0,N,nbf)

initF <- matrix(0,T,nbf)
for (k in 1:nbf) {
initF[k,k] <- sqrt(T)
}

beta <- coefsbai(var$y,var$X,nbf,initL,initF,iter,1)
sbeta<- stdbai_iid(var$y,var$X,beta$coefX,beta$coefC,beta$coefi,
                   beta$coeft,beta$coefL,beta$coefF)

expl <- 1-(ncol(var$X)==1)
if (expl==0) { 
	cbai[s,] <- beta$coefX 
	vbai[s,] <- sbeta^2
}
if (expl==1) { 
	cbai[s,1] <- t(beta$coefX)[1,1] 
	vbai[s,1] <- sbeta[1,1]^2 
}
}

# SORTIE D'UN TABLEAU DE RESULTATS

results <- matrix(0,3,1)

# R�sultats Bai (2009)

results[1,1] <- mean(cbai[,1])
results[2,1] <- var(cbai[,1])
results[3,1] <- mean(vbai[,1])

results
}





##########################################################################
# SIMULATIONS 

N 		<- 143
N0 		<- 13
T		<- 20
T0		<- 8

S <- 1000

nl <- c("Interactive effects,","linear counterfactual","","Interactive effects,","linear","",
	  "Interactive effects,","matching","","Interactive effects,","constrained","",
        "Synthetic controls","","","First difference","","")

iter <- 20
iterEM <- 1

dfei0  <- read.table(file="C:/fei0.txt")
dfet   <- read.table(file="C:/fet.txt")
dvecL0 <- read.table(file="C:/vecL0.txt")
dvecF  <- read.table(file="C:/vecF.txt")

fei0	 <- as.matrix(c(dfei0[1:N0,1],dfei0[(nrow(dfei0)-N+N0+1):nrow(dfei0),1]),N,1)
fet    <- as.matrix(dfet[1:T,],T,1)

vecL0	 <- as.matrix(rbind(dvecL0[1:N0,],dvecL0[(nrow(dvecL0)-N+N0+1):nrow(dvecL0),]),N,ncol(dvecL0))
vecF	 <- as.matrix(dvecF[1:T,],T,nbf)



# TABLE 1, PAPER # 
##################

link <- "normal"
nbf <- 2

nc <- c("0",".5","1")
lbl <- "Support difference"

res <- matrix(0,18,3)

supp <- 0
results <- simul2(N,N0,T,T0,alphaT,alphaX,alphaC,supp,
                  fei0,fet,vecL0,vecF,iter,iterEM,link,nbf,S,0,0,1)
res[,1] <- results 

supp <- .5
results <- simul2(N,N0,T,T0,alphaT,alphaX,alphaC,supp,
                  fei0,fet,vecL0,vecF,iter,iterEM,link,nbf,S,0,0,1)
res[,2] <- results 

supp <- 1
results <- simul2(N,N0,T,T0,alphaT,alphaX,alphaC,supp,
                  fei0,fet,vecL0,vecF,iter,iterEM,link,nbf,S,0,0,1)
res[,3] <- results 

textable(x=res,
        file="C:\\table1_paper",
        label=lbl,nomligne=nl,nomcol=nc, append=FALSE,decim=3,) 



# TABLE 2, PAPER - sinusoidal trend #
#####################################

# First factor becomes a sinusoidal trend

pi=3.14159265359

vecF	 <- as.matrix(dvecF[1:T,],T,nbf)
vecF[1:T,1] <- 5*sin(pi*as.matrix(seq(1,T),T,1)/T)

link <- "normal"
nbf <- 2

nc <- c("0",".5","1")
lbl <- "Support difference"

res <- matrix(0,18,3)

supp <- 0
results <- simul2(N,N0,T,T0,alphaT,alphaX,alphaC,supp,
                  fei0,fet,vecL0,vecF,iter,iterEM,link,nbf,S,0,0,1)
res[,1] <- results 

supp <- .5
results <- simul2(N,N0,T,T0,alphaT,alphaX,alphaC,supp,
                  fei0,fet,vecL0,vecF,iter,iterEM,link,nbf,S,0,0,1)
res[,2] <- results 

supp <- 1
results <- simul2(N,N0,T,T0,alphaT,alphaX,alphaC,supp,
                  fei0,fet,vecL0,vecF,iter,iterEM,link,nbf,S,0,0,1)
res[,3] <- results 

textable(x=res,
        file="C:\\table2_paper",
        label=lbl,nomligne=nl,nomcol=nc, append=FALSE,decim=3,) 



# TABLE 3, PAPER - number of factors # 
######################################

vecF	 <- as.matrix(dvecF[1:T,],T,nbf)

link <- "normal"
supp <- 0

nc  <- c("1","2","3","4","5")
lbl <- "Number of factors"

res <- matrix(0,18,5)

for (nbf in 1:5) {

results <- simul2(N,N0,T,T0,alphaT,alphaX,alphaC,supp,
                  fei0,fet,vecL0,vecF,iter,iterEM,link,nbf,S,0,0,1)

res[,nbf] <- results 
}

textable(x=res,
        file="C:\\table3_paper",
        label=lbl,nomligne=nl,nomcol=nc, append=FALSE,decim=3,) 



# TABLE B1, PAPER - uniform errors # 
####################################

link <- "uniform"
nbf <- 2

nc <- c("0",".5","1")
lbl <- "Support difference"

res <- matrix(0,18,3)

supp <- 0
results <- simul2(N,N0,T,T0,alphaT,alphaX,alphaC,supp,
                  fei0,fet,vecL0,vecF,iter,iterEM,link,nbf,S,0,0,1)
res[,1] <- results 

supp <- .5
results <- simul2(N,N0,T,T0,alphaT,alphaX,alphaC,supp,
                  fei0,fet,vecL0,vecF,iter,iterEM,link,nbf,S,0,0,1)
res[,2] <- results 

supp <- 1
results <- simul2(N,N0,T,T0,alphaT,alphaX,alphaC,supp,
                  fei0,fet,vecL0,vecF,iter,iterEM,link,nbf,S,0,0,1)
res[,3] <- results 

textable(x=res,
        file="C:\\tableB1_paper",
        label=lbl,nomligne=nl,nomcol=nc, append=FALSE,decim=3,) 





# TABLE B2, PAPER - number of periods # 
#######################################

link <- "normal"
supp <- 0
nbf <- 2

nc <- c("T=20,T0=8","T=10,T0=4")
lbl <- "Number of periods"

res <- matrix(0,18,2)

results <- simul2(N,N0,T,T0,alphaT,alphaX,alphaC,supp,
                  fei0,fet,vecL0,vecF,iter,iterEM,link,nbf,S,0,0,1)

res[,1] <- results 

fet2    <- as.matrix(dfet[5:14,],10,1)
vecF2	  <- as.matrix(dvecF[5:14,],10,nbf)

T		<- 10
T0		<- 4

results <- simul2(N,N0,T,T0,alphaT,alphaX,alphaC,supp,
                  fei0,fet2,vecL0,vecF2,iter,iterEM,link,nbf,S,0,0,1)

T		<- 20
T0		<- 8

res[,2] <- results 

textable(x=res,
        file="C:\\tableB2_paper",
        label=lbl,nomligne=nl,nomcol=nc, append=FALSE,decim=3,) 



# TABLE B3, number of units # 
#############################

link <- "normal"
supp <- 0
nbf <- 2

nc <- c("N=143,N0=13","N=275,N0=25")
lbl <- "Number of individuals"

res <- matrix(0,18,2)

results <- simul2(N,N0,T,T0,alphaT,alphaX,alphaC,supp,
                  fei0,fet,vecL0,vecF,iter,iterEM,link,nbf,S,0,0,1)

res[,1] <- results 

N 		<- 275
N0 		<- 25

fei02	 <- as.matrix(c(dfei0[1:N0,1],dfei0[(nrow(dfei0)-N+N0+1):nrow(dfei0),1]),N,1)
vecL02 <- as.matrix(rbind(dvecL0[1:N0,],dvecL0[(nrow(dvecL0)-N+N0+1):nrow(dvecL0),]),N,ncol(dvecL0))

results <- simul2(N,N0,T,T0,alphaT,alphaX,alphaC,supp,
                  fei02,fet,vecL02,vecF,iter,iterEM,link,nbf,S,0,0,1)

N 		<- 143
N0 		<- 13

res[,2] <- results 

textable(x=res,
        file="C:\\tableB3_paper",
        label=lbl,nomligne=nl,nomcol=nc, append=FALSE,decim=3,) 



# TABLE B4, SERIAL TIME AUTOCORRELATION # 
#########################################

link <- "normal"
supp <- 0
nbf  <- 2

nc <- c(".1",".3",".5",".7",".9")
lbl <- "Support difference"

res <- matrix(0,18,5)

rho  <- 0.1
results <- simul2(N,N0,T,T0,alphaT,alphaX,alphaC,supp,
                  fei0,fet,vecL0,vecF,iter,iterEM,link,nbf,S,rho,0,2)
res[,1] <- results 

rho  <- 0.3
results <- simul2(N,N0,T,T0,alphaT,alphaX,alphaC,supp,
                  fei0,fet,vecL0,vecF,iter,iterEM,link,nbf,S,rho,0,2)
res[,2] <- results 

rho  <- 0.5
results <- simul2(N,N0,T,T0,alphaT,alphaX,alphaC,supp,
                  fei0,fet,vecL0,vecF,iter,iterEM,link,nbf,S,rho,0,2)
res[,3] <- results 

rho  <- 0.7
results <- simul2(N,N0,T,T0,alphaT,alphaX,alphaC,supp,
                  fei0,fet,vecL0,vecF,iter,iterEM,link,nbf,S,rho,0,2)
res[,4] <- results 

rho  <- 0.9
results <- simul2(N,N0,T,T0,alphaT,alphaX,alphaC,supp,
                  fei0,fet,vecL0,vecF,iter,iterEM,link,nbf,S,rho,0,2)
res[,5] <- results 


textable(x=res,
        file="C:\\tableB4_paper",
        label=lbl,nomligne=nl,nomcol=nc, append=FALSE,decim=3,) 





# TABLE B5, HETEROGENEOUS VARIANCES # 
#####################################

link <- "normal"
supp <- 0
nbf  <- 2

nc <- c("1.2","1.5","2","2.5","3")
lbl <- "Variance ratio"

res <- matrix(0,18,5)

Vr <- 1.2
results <- simul2(N,N0,T,T0,alphaT,alphaX,alphaC,supp,
                  fei0,fet,vecL0,vecF,iter,iterEM,link,nbf,S,0,Vr,3)
res[,1] <- results 

Vr <- 1.5
results <- simul2(N,N0,T,T0,alphaT,alphaX,alphaC,supp,
                  fei0,fet,vecL0,vecF,iter,iterEM,link,nbf,S,0,Vr,3)
res[,2] <- results 

Vr <- 2
results <- simul2(N,N0,T,T0,alphaT,alphaX,alphaC,supp,
                  fei0,fet,vecL0,vecF,iter,iterEM,link,nbf,S,0,Vr,3)
res[,3] <- results 

Vr <- 2.5
results <- simul2(N,N0,T,T0,alphaT,alphaX,alphaC,supp,
                  fei0,fet,vecL0,vecF,iter,iterEM,link,nbf,S,0,Vr,3)
res[,4] <- results 

Vr <- 3
results <- simul2(N,N0,T,T0,alphaT,alphaX,alphaC,supp,
                  fei0,fet,vecL0,vecF,iter,iterEM,link,nbf,S,0,Vr,3)
res[,5] <- results 

textable(x=res,
        file="C:\\tableB5_paper",
        label=lbl,nomligne=nl,nomcol=nc, append=FALSE,decim=3,) 



# TABLE B6, FINITE SAMPLE AND ASYMPTOTIC VARIANCES FOR BAI (2009) # 
###################################################################

nl <- c("Mean coeff","Std coeff","Square root Asymp. Var")

link <- "normal"
supp <- 0

nc  <- c("1","2","3","4","5")
lbl <- "Number of factors"
 
res <- matrix(0,3,5)

for (nbf in 1:5) {

results <- simul3(N,N0,T,T0,alphaT,alphaX,alphaC,supp,
                  fei0,fet,vecL0,vecF,iter,iterEM,link,nbf,S)

results[2,] <- sqrt(results[2,])
results[3,] <- sqrt(results[3,])

res[,nbf] <- results 
}

res[1,] <- res[1,]-.3
res[1,] <- round(res[1,]*1000)/1000

textable(x=res,
        file="C:\\tableB6_paper",
        label=lbl,nomligne=nl,nomcol=nc, append=FALSE,decim=3,) 




